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The motion of an ion in a coherent lower hybrid wave (characterized by |fc||| <C 
and w 3> rii) in a tokamak plasma is studied. For ions satisfying v± > uj/k^, the 
Lorentz force law for the ions is reduced to a set of difference equations which give the 
Larmor radius and phase of an ion on one cyclotron orbit in terms of these quantities a 
cyclotron period earlier. From these difference equations an earlier result [Phys. Fluids 
21, 1584 (1978)] that above a certain wave amplitude the ion motion is stochastic, is 
readily obtained. The stochasticity threshold is given a simple physical interpretation. 
In addition, the difference equations are used to derive a diffusion equation governing 
the heating of the ions above the stochasticity threshold. By including the effects of 
collisions, the heating rate for the bulk ions is obtained. 
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I. INTRODUCTION 

The motion of an ion in a plasma in the presence of a lower hybrid wave becomes 
stochastic if the amplitude of tlie lower hybrid wave exceeds a thrcsliold. Tliis provides 
a mechanism by which ions may be directly heated by the wave when using lower hybrid 
waves to heat a tokamak plasma. In this paper we extend the work of an earlier papor^ 
(henceforth referred to as I), in which the stochasticity threshold was derived. The main 
objects are to give a simple physical explanation of the stochasticity threshold, to derive the 
velocity space diffusion coefficient, and hence to determine the heating rate. 

In I it was shown that for the purposes of computing the ion motion it may be assumed 
that the magnetic field is uniform if the fractional change in the magnetic field over a 
Larmor radius is small. Similarly, the lower hybrid wave may be treated as perpendicularly 
propagating if the wavevector, k, satisfies |A;[|/fcx| < j^i)^^!^ . Both these conditions 
are well satisfied in normal circumstances. Shear in the magnetic field may also affect the 
motion;^ however, if the scale length of the shear greatly exceeds the Larmor radius, we 
expect that this effect may also be ignored. Thus, in this paper, we consider only the case 
of an electrostatic wave propagating perpendicularly to a uniform magnetic field; i.e., 

B = Boz, E = £;oycos(fc±y -t^t). (1) 

As in I these fields are taken to be imposed, since only the motion of a small fraction of the 
ions on the tail of the ion distribution function becomes stochastic, while the bulk ions and 
electrons support the wave. Having determined the heating rate it will be possible to allow 
for the damping of the wave. 

Following I the Lorentz force law for an ion with mass, m^, and charge, qi — ZiC, may 
be written as 

y + y = acos{y-vt), x = y, (2) 
where length is normalized to fcj^ and time to (fij = qiBo/nii) and where 

u = co/Qi, (3a) 

Eq/Bq 

a = — . 

For later use we shall also define a normalized Larmor radius, r, by 

r = k±v±/fli, (3c) 

where v"^ = ir? + y"^ . 

The plan of this paper is as follows: We first reduce the Lorentz force law for the ion, 
(2), to a set of difference equations (Sec. II) which determine the velocity and phase of the 
ion in terms of its velocity and phase a cyclotron period earlier. These difference equations 
have several symmetries which are presented in Sec. Ill and which simplify the study of the 
ion motion. This is followed by a study of the ion motion for infinitesimal i?o (Sec. IV) 
where we show that the motion when w is at a harmonic or half harmonic of the cyclotron 
frequency has a very different character from the motion at other frequencies. In Sec. V 
we describe qualitatively the motion when Eq is finite. Here we define the stochasticity 
threshold and explain it physically. The statisical properties of the ion motion above the 
stochasticity threshold are described by the Krylov-Kolmogorov- Sinai entropy and the cor- 
relation function which are examined in Sees. VI and VII. From the correlation function we 
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find the diffusion coc;f5ci«it (Sc;c;. VIII). For times much longer than the correlation time 
the ion motion in perpendicular velocity space is governed by a diffusion equation which 
is derived in Sec. IX. This equation is checked against the exact equations of motion by 
comparing a Monte Carlo solution of the diffusion equation with a simulation described in 
I. Since the stochastic heating only affects the tail ions we include a collision term in the 
diffusion equation as a means of heating the bulk ions and electrons. The resulting Fokker- 
Planck equation which describes the ion motion in two-dimensional velocity space {v^_ and 
v\\) is reduced to a one-dimensional equation in va_ only (Sec. X). The one-dimensional Fok- 
ker-Planck equation is checked against the two-dimensional equation by numerically solving 
each of them (Sec. XI). Analytic steady-state solutions to the one-dimensional equation are 
derived in Sec. XII which give results for the heating rates of the ions and electrons. 



II. DERIVATION OF THE DIFFERENCE EQUATIONS 

In this section we derive difference equations which approximately describe the motion 
of the ion. The procedure involves integrating the equations of motion along unperturbed 
orbits. Nonlinearity is introduced because we repeatedly correct the unperturbed trajecto- 
ries. We assume that the wave frequency is much larger than the cyclotron frequency; i.e., 
ljj ':$> Vli or V \ . (This is the case for lower hybrid waves.) On a time scale between the 
wave period and the cyclotron period the ion does not experience the effects of the magnetic 
field and behaves, in some respects, as though in zero magnetic field. Thus, when it passes 
through the wave-particle resonance points, o) = kv or y = v (Fig. 1), it will experience 
kicks, which we can approximate by impulses. This is the physical picture that will guide 
our derivation. We shall not be too concerned yet about the limits in which our approx- 
imations are valid, prefering to rely on physical intuition. Appendix A gives an alternate 
derivation starting from the Hamiltonian. While this derivation is less transparent, it does 
enable us to state the limits of validity of the difference equations. 

From Fig. 1 we sec that for r > v the particle passes through wave-particle resonance 
twice per cyclotron orbit. We will approximate the force on the particle due to the electric 
field by impulses or delta functions at these two points and zero elsewhere. 

In order to evaluate the area of the delta function forces we Taylor-expand the orbit 
about the wave-particle resonance point using the unperturbed equation of motion, y+y = 0. 
Thus, we have 



where t' = t — tc, 4>c = Vc — i^tc, tc and yc are the time and position of the wave-particle 
"collision," and, by assumption, y{tc) = v. We substitute (4) into the electrostatic force 
term in (2), a cos(y — vt)^ and approximate the resulting expression by B 5{t — tc), where 6 
is the Dirac delta function and 



y 



i't = yc + y{tc)t' + lyitc)^^ - vt 



= 0c - ^yc^'^, 



(4) 



B = 




Q:(27I"/ Ij/cI)^^^ cos[0c - sign(2/c)7r/4] 



(5) 
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[sec Rcf. 3, Eq. (4.3.144)]. Wc note that (6,; and hence B are the same (to first order in 
a) whether we compute tc and yc from the intersection of the orbit for t < tc with y = v 
or from the intersection of the orbit for t > tc with y = u. This is so because yc and tc 
are unchanged by the kick the ion receives, and because d{y — vt)/dt = near the collision 
(since y = v). As a consequence the mapping we define with the difference equations will 
be invariant to time reversal, which in turn means that the mapping is measure-preserving, 
a necessary condition for a non-dissipitive (i.e., Hamiltonian) system. 

We define the cyclotron orbit as beginning and ending at y = 0, ?) < (Fig. 1). 
Subscript, j, will refer to the orbit at the beginning of the jth or ending of the [j — l)th 
orbit. A subscript, j + i, will be used to refer to the middle of the jth orbit, where y = 0, 
y > 0. Finally, subscripts, + and — , will refer to the collisions with j/c > and yc < 0. 

We describe the orbit by two parameters, a reduced Larmor radius, 

p={r'^- v^yl'^ - vcos~'^{v/r) + - 7r/4, (6a) 
and the phase of the wave, 

9 = iyt (mod 27r). (6b) 

(Note that tor r ^ u, p r.) It is convenient to introduce the sum and difference of these 
quantities, v = 6 + p, u = 9 — p. 

Evaluating B_ we note that (f>c = —{fj— v^Y^"^ ~ ^fe + ^ ~ cos~^(i//rj)] = —Vj — 7r/4, 
so that 

B_ = a(27r)i/2(r2 _ i.2)-i/4 ^^g^^.^ 

In (7) we have intentionally dropped the subscripts on r in the amplitude term. The 
justification for this is that r is large, so a change of order a in r has negligible effect on the 
amplitude factor compared with the same change in the argument to the cosine. 
Working to order a we then find 

rj+i/2^rj + {u/r)B_, (8a) 

tj+1/2 =tj+n- ^ , (8b) 

whence from (6) 

Pj+i/2 = Pj + ^{r'' - uy/\-^B_, (9a) 

^^.+1/2 = Oj + i/TT - u{r^ - u'^f'^r-^B_ (9b) 
[note that dp/dr — (r^ — i/^)!/^/^], and from (5) 

/2 = ■Wj + !^7r — 21^ A cos Vj , (10a) 
Vj+l/2 = Vj + VTV, (10b) 



where 



^ ,2y'^»«f-^^ (11) 
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Similarly, we may show that 

Uj+l = Uj+l/2 + z^TT, (12a) 

Vj+i = + UTT + 27r^cosWj+i. (12b) 

Combining (10) and (12) we obtain the following difference equations which give the (j+l)th 
orbit in terms of the jth orbit, 

u = 6 — p, V = + p, (13a) 

e=^{v + u), p = ^{v-u), (13b) 

Uj+i — Uj = 2^6 — 27rAcosvj, (13c) 

Vj+i — Vj = 2tt6 + 27rAcosMj+i, (13d) 



where 

av 
A= — 
r 



ff(i)'(r) , (14) 
6 = u — n, (15) 

H^}^ is the Hankel function of the first kind and of order \S\ < ^, and n is an integer. In 
writing (14) we have made use of the expansion of the Bcsscl functions for r > z/ + {^I'Y^^ 
[Rcf. 3, Eq. (9.3.3)]. Writing A in this way is mainly a notational convenience, although it 
comes out naturally in this form when we use the whole cyclotron orbit as the unperturbed 
orbit.* Wc; treat A as being a parameter (independent of p) in (13). We shall define T as 
the mapping which takes a point to its iterate: 

Tie„p,) = {9,+uPj+i). (16) 
The limits of validity of (13) are found to be (see Appendix A) 

i/»l, r-^/»(i^y)V^ A ^ {r'' - u^/^r^ (17) 

Note that the problem now just depends on two parameters, 6 and A. This was rec- 
ognized in I, but was only established analytically in the case r ^ u, which is much more 
restrictive than (17). 

The difference equations, (13), may be iterated numerically and compared with the 
solution to the Lorentz force law, (2). Such a comparison is made in Fig. 2, where Fig. 
2(a) is taken from 1 and Fig. 2(b) is obtained from the difference equations. We see that 
there is excellent agreement, so that the difference equations indeed provide an accurate 
approximation to the Lorentz force law. 

We note that the kicks received by the ions happen at the Landau resonance points 
(Fig. 1), and so it is interesting to compare the ion motion under the conditions described 
by (17) with Landau damping in the absence of a magnetic field. There are two important 
differences. Because the magnetic field sweeps the vector, v_l, through all angles, the res- 
onant particles in our case are those for which r > u ov v± > ui/k±. These particles are 
much more numerous than those satisfying the normal Landau resonance condition y = v ov 
Vy = Lo/k±. The second difference is that the kick that the particle receives, B, is a function 
of the magnetic field, since the magnetic field causes the particle to spend only a short time 
in resonance. In the linear limit of normal Landau damping the resonant particles remain 
in resonance forever. 
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III. SIMPLE PROPERTIES OF THE DIFFERENCE EQUATIONS 

The difFcrcncc equations, (13), describe the mapping of the {0,p) plane onto itself. It 
models the nonlinear coupling of two harmonic oscillators, the parameter, (5, giving the 
relation between the unperturbed frequencies and A giving the strength of the coupling. 
The mapping is area-preserving, being derivable from the generating function, 

Vj) = Uj+iVj + 2'K5(uj+i — Vj) + 27rA(sin?ij+i + sinwj). (18) 

[Equations (13c) and (13d) are given by uj = dF/dvj and Vj+i = dF / duj+i, respectively.] 
Thus, (13) describes a Hamiltonian system. 

The [9, p) plane is periodic (period 27r) in both the 9 and p directions and so is topo- 
logically equivalent to a 27r x 27r torus. The periodicity in p is a consequence of taking A to 
be independent of p. This, in turn, will allow a very simple determination of the diffusion 
coefficient (see Sees. VII and VIII). The periodicity in u and v allows a stronger statement 
to be made, viz., the transformation, 

9^9 + t:, p->p±7r, (19) 

leaves (13) invariant. This means it is sufficient to examine a range of 27r in 9 and tt in p. 
Other transformations that leave the difference equations invariant are: 

5^5 + 1; (20) 

5^-5, 9 ^TT-e, p -^-p; (21) 

A^-A, e ^9--k; (22) 

3 ^-j, 6 ^TT- 9. (23) 

The first three of these allow a further restriction of the problem to < (5 < | and A > 0. 
The last transformation shows the invariance of (13) to time reversal. 



IV. SMALL AMPLITUDE SOLUTION 

In this and the next section wc examine some aspects of the transition from coherent 
to stochastic behavior, beginning with the analysis for small A. 

It is possible to write down the Hamiltonian for the difference equations, (A13). For 
small A, canonical transformations similar to those in I may be performed, which allow a 
solution to the problem. Here we take a different approach which achieves the same results 
in a more direct manner. 

In the limit, ^ = 0, (13) describes two uncoupled harmonic oscillators. If (5 is a rational 
number, s/p, then all the points in the {9,p) plane are pth-order fixed points. If we now 
let A be small, we expect 9p — 9q and pp — pq to be small also. In that case the difference 
equations may be approximated be differential equations for 9 and p. (It might seem that we 
are going round in circles here. However, these differential equations will be much simpler 
than the exact equations of motion, since they describe motion with only one degree of 
freedom.) If we can find the Hamiltonian, h, such that 9 = dh/dp and p = —dh/d9, then 
we have found an additional constant of the motion, h. Lines of constant h{9, p) give the 
trajectories of ions in the {9, p) plane. 
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We begin with the case, 6 = (i.e., a cyclotron harmonic). Then (13) becomes 

wi — Wo = — 27rAcost;o, (24a) 

vi-vo = 27rAcosui = 27r>lcosuo + 0{A'^), (24b) 

so that 

01-00 =2'KAsmeosmpo + 0{A^), (25a) 

pi- Po = 2'kAcos0ocospo + O{A^). (25b) 
These equations are Euler's approximation to the solution of the differential equations, 

= AsmOsmp, (26a) 

p = A cos cos p. (26b) 

[We define 0j = 0{t = 27rj), etc.] These, in turn, are Ifamilton's equations for the Hamilto- 
nian, 

h = —Acospsm0. (27) 

The trajectories of the ion in {0,p) space, given by this Hamiltonian, are shown in Fig. 3(a). 
This result is the same as that obtained by Timofeev^ using the full Hamiltonian for the 
ion. Note that all particles exchange energy with the wave, but that the amount of energy 
exchanged is bounded in time. 

It is instructive to compare this result with standard linear theory. As we noted, our 
approach is nonlinear because we repeatedly correct the unperturbed orbits. If this is not 
done (almost) all the orbits are secular with pjy — po ^ AN . These orbits would appear as 
straight lines in the {0,p) plane which are tangent to the curves drawn in Fig. 3(a) at the 
initial positions of the ions. Because all the orbits arc secular, the linear dielectric function 
is divergent at a cyclotron harmonic (for fcy = 0). Linear theory is valid for t < t where t is 
the inverse frequency of motion around the islands in Fig. 3(a). From (27) we find r = A~^ 
(near the center of the islands). 

For 6 = ^ (i.e., a; at a half harmonic of Oj) we find 

U2-uo = 2tt + 47r^A^ sin vo cos wo + 0{A^), (28a) 

W2 - i;o = 27r - Att'^A'^ coswo sinwo + 0{A^). (28b) 
Ignoring the shift of u and v by 27r the differential equations for and p are 

= iTA^sm2p, (29a) 

p= -7r^^sin26», (29b) 

which are derivable from the Hamiltonian, 

h = -nA'^cos{0 + p)cos{e-p). (30) 

The trajectories for this case are shown in Fig. 3(b). 

The motion is very similar to the case, 6 = 0, Fig. 3(a). The most important difference 
is the extra factor of A in (30) as compared with (27) so that the exchange of energy with 
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thc wave is much slower at a half harmonic. Indeed, to order A there is no exchange of 
energy since the phase of the ion relative to the wave increases by about tt each cyclotron 
period, and so the kicks received by the ion on successive cyclotron orbits nearly cancel. 
Thus, there is no contribution to the linear dielectric function from the half harmonics. We 
expect the 0{A^) change in the speeds of the particles in this case to lead to additional 
nonlinear terms in the dielectric function which would cause a nonlinear damping of the 
waves due to half-harmonic resonances. 

In the general case where 5 = s/p ^Q,^,we have 

0p-0o=27rs+p7r2A2^^^%-^+O(A3), (31a) 

sm(7ro) 

Pp-po = 0{A^). (31b) 

These expressions were explicitly checked only for 5 = g, i, |, The algebra 
required to do the iterations and expansions of the difference equations is quite onerous, and 
so was performed using the algebraic manipulation system, MACSYMA.^ We note that the 
expression for 6 given by dividing (31a) by 2'Kp agrees with the expression for the frequencies 
given by Eq. (43) of I. [This result in I was derived using canonical transformations on the 
the Hamiltonian, (Al), and required considerably more effort than (31).] The trajectories in 
this case are p = const. [Fig. 3(c)]. In this case there is no net exchange of energy between 
the particles and the wave. 



V. THE STOCHASTIC TRANSITION 

For finite A the motion becomes much more complicated. Here we briefly discuss what 
happens. The exposition is intended to complement that given in I, although there is some 
overlap. General reviews of the transition to stochasticity may be found in Refs. 7 9. 

We begin by noting the topological difference between the cases 6 = 0,^ and 5^0, 
5 (Fig. 3). This leads to differences in the way that stochasticity develops, and because of 
this we define the stochasticity threshold differently in the two cases. 

With S not close to or i, higher-order terms in the difference equations, i.e., the 
0{A^) terms in (31), cause islands to appear. The location of the islands may be derived by 
requiring that (9p — 6'o)/(27rp), (31a), be rational. Where these islands overlap the motion 
becomes stochastic [Fig. 2(b)]. However, the KAM (Kolmogorov-Arnold-Moser) theorem^" 
ensures that for sufficiently small A some trajectories exist which span 6 [again see Fig. 
2(b)]. [These are KAM "surfaces" reflecting their dimensionality in the original problem, 
although they only appear as lines in the {9, p) plane.] These act as barriers, preventing p for 
a given particle from increasing or decreasing by more that about tt (since the KAM surfaces 
repeat periodically in p). At some value of A the last KAM surface disappears, allowing 
unrestriced motion in p. We deflne this value of A to be the stochasticity threshold, Ag. In 
I the value of Ag was found to be about |. Recently accurate methods for determining As 
have been developed by Greene, -'^'^ although applying these methods to the problem in hand 
presents difficulties because the relevant resonances do not exist down to A = 0. 

For 5 equal to or ^, we have, for small A, all of phase space covered by islands [Fig. 
3(a) and (b)]. When S is not equal to or i but is close to these values, Fig. 3(c) applies 
only for very small amplitudes. If ^ > \S\ or A > (|1 — 2S\ /tt)^/^, first- and second-order 
islands appear. At larger amplitudes these islands can cover nearly all of phase space giving 
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a picture similar to Fig. 3(a) and (b). The effect of finite A in these cases is to prockice 
chains of islands within the main islands. These overlap close to the separatrix causing a 
layer around the separatrix to become stochastic. For 5 = the thickness of this layer is 
found^^ to be on the order of 

exp[-7r/(2A)] (32) 

for small A\ a similar exponential dependence is expected for 5 = ^. Thus, in contrast with 
the previous case it is possible for ions to be heated even at low amplitudes. However, for 
small A, the thickness of the stochastic layer, (32), is extremely small. This means firstly 
that only a tiny fraction of particles will be in the stochastic region and secondly that those 
that are will be heated very slowly. Therefore we define the stochasticity threshold in this 
case as the value of A, Ag, for which the stochastic layer occupies a substantial fraction of 
phase space. Note that this docs not define a precise value of As- The definition is useful 
however since (32) is a strong function of A. In I it was found that As had approximately 
the same value as for 5 7^ 0, \, namely \. 

The stochasticity threshold. As = \, has a simple physical interpretation. At this 
value of A, the kick received by the particle on one transit through wave-particle resonance 
(Fig. 1) is sufficient to change the phase that the particle sees when next in wave-particle 
resonance by 7r/2 ( = 211 As). This explanation, which is valid for r > v, complements the 
explanation in terms of trapping for the case r k, v, which was given in I. 

We end this section by looking at the eigencurves of the mapping. These are curves 
in the (^, p) plane which map into themselves on applying the mapping, T . A study of the 
eigencurves provides another view of the stochastic transition. Again the reader is referred 
to Refs. 7-9 for a fuller discussion of this. We consider only the case i5 = 0. (Because of the 
additional symmetries in this case, the eigencurves are easier to generate. Similar behavior 
is observed for other values of 5.) When A is infinitesimal the eigencurves are given by 
Fig. 3(a); and Fig. 4 shows them for finite values of A. For finite but small A, most of 
the eigencurves remain closed; these may be generated by iterating the mapping of a single 
point many times. However, the eigencurves emanating from the hyperbolic fixed points no 
longer meet up. They are now infinitely long open curves. They are generated as follows: 
Close to the hyperbolic fixed point, the mapping may be linearized. The eigencurves of the 
linearized mapping are hyperbolae. The asymptotes of the hyperbolae are eigencurves on 
which a point either moves away from or towards the fixed point. These eigencurves are 
called the unstable and stable manifolds respectively. Short straight line segments, which 
connect a point with its image under T . are chosen on these ^instable and stable manifolds. 
These line segments arc then repeatedly mapped forwards and backwards by applying T and 
T~^; and the union of the mapped segments generates the open eigencurves. In practice, 
only a few (less than 10) iterations are required. Figure 4(a) shows the eigencurves for 
A = 0.3. We see that the open curves intersect each other at a finite angle, and start 
oscillating as they approach the fixed point. Because the mapping is area-preserving these 
oscillations grow larger and larger as the fixed point is approached. These curves occupy a 
finite area which may be shown to be stochastic. As illustration of this we show the iteration 
of a single point starting close to the hyperbolic fixed point (Fig. 5). We see that it covers 
nearly all the area outside the first-order islands. All these points lie on a single eigencurve. 
The area occupied by the eigencurves (i.e., the area of phase space that is stochastic) is 
proportional to the angle at which the eigencurves first intersect (approximately half way 
between the; fixed points) . This then gives another way, in addition to the method of island 
overlap, of estimating the fraction of phase space that is stochastic. This angle again has 
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thc exponential dependence on amplitude given in (32). In Fig. 4(b), A has been increased 
to 0.35. Here we see the open eigencurves intersecting at a larger angle. In addition, the 
elliptic fixed points in Fig. 4(a) have turned into hyperbolic fixed points with reflection [from 
(B7) we find that this happened aX, A = tt~^ w 0.318]. Thus, open eigencurves emanate 
from these points. However, the angle at which they intersect is so small that the area they 
occupy is comparable to the thickness of the lines used to draw the figure. Eigencurves 
corresponding to the original islands are still present. However, these have disappeared at 
A = 0.45, Fig. 4(c). At this amplitude the two sets of open eigencurves intersect each other. 
This may be seen by comparing this figure with Fig. 4(d), where we have extended one of 
the eigencurves emanating from the hyperbolic fixed point with refiection. 



VI. THE KRYLOV-KOLMOGOROV-SINAI ENTROPY 

Having established that the ion motion becomes stochastic for A > Ag, we need to be 
able to describe the motion above the stochasticity threshold. One of the most important 
parameters in this regard is the Krylov-Kolmogorov-Sinai (KS) entropy, h. This is a mea- 
sure of the local instability of trajectories and is defined as the average rate of divergence 
of neighboring (infinitely close) trajectories. Thus, after N iterations the distance between 
neighboring particles initially do apart is approximately dQexp{hN) (for do small and N 
large). 

The KS entropy is important for two reasons. Firstly, given a group of particles which 
initially occupy a small region of phase space, Ap, we can estimate the number of itera- 
tions for the phases, 6, of the particles to become random to be 

log(27r/A6i)//i. (33) 

Secondly, the exponential divergence of trajectories leads to a mixing of phase space, which 
in turn results in the decay of the correlation function, Ck (Sec VII), and allows the motion 
of the particle to be described by a diffusion equation. 

In order to determine h we define a reference trajectory, {uj, Vj). [It is more convenient 
to work in {u,v) space.] We consider a neighboring trajectory, {uj + Suj,Vj + 5vj), where 
6uj and 6vj are infinitesimals. Then 5uj and 5vj satisfy 

where 

/ 1 27rAsini),- \ 

■•^•^ o . • 1 . 2,2 . . • (35) 

K—ZTrAsmUj+i 1 — 47r A smuj+ismvj J 



If we define M jv by 



JV-l 



then h is given by 



3=0 

h= Lim(log|A;v|/A^), (37) 

Af— >oo 

where Ajv is the largest (in magnitude) eigenvalue of Mjv- This definition is equivalent to^ 

h = (log(£,+i/^,)), (38) 
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where £j is the length of the vector {5uj, 5vj) and the average is taken over a particle orbit. 

Figure 6 shows log|AAr|/A^ as a function of n for various values of A. We see that 
it does have a limiting value as ^ oc' and so h is well defined. Figure 7 shows h as 
a function of ^ for i5 = 0.23. For A small, h is close to zero, because when the motion is 
coherent neighboring trajectories diverge linc^arly. rather than exponentially. When A is still 
below As, h becomes finite, since the reference trajectory was chosen in the stochastic part 
of phase space (Fig. 2). As A passes through the stochasticity threshold there is a slight 
drop in h. Finally, for A ^ Ag, h is on the order of log A. 

We may derive the expression for h in the limit A^ Ag. If and Hj are the larger 
and smaller (in magnitude) eigenvalues of J j , then for A^ Ag we have 

Aj « — 47r^A^ sinuj+i sinwj, ^j = l/Xj. (39) 

Furthermore, the eigenvectors corresponding to these eigenvalues are approximately v and 
u respectively. Since the eigenvectors for different Jj are nearly parallel to one another, the 
eigenvalues of the product of the Jj are approximately the product of the eigenvalues of the 
Jj] i.e., 

N-l 

An^UXj. (40) 

Substituting (39) into (40) we have 

h = (log|Aj|) = 21og(7rA) + (log|4sinUj-+isinz;j|), (41) 

where the average is taken over a particle orbit. Finally, we note that for ^ » Ag, a 
particle's orbit wanders over most of phase space, spending equal amounts of time in equal 
areas of phase space (i.e., the particle's orbit is approximately ergodic). Thus, the average 
along the orbit may be replaced by a phase space average. The last term in (41) becomes 

2(27r)-W log|2sin^| # = (42) 
Jo 

[see Ref. 3, Eq. (4.3.145)]. Thus, for A > A^, we have 

/i w 21og(7rA). (43) 
This is shown as a dashed line in Fig. 7. 



VII. THE CORRELATION FUNCTION 

We saw in the last sections that a group of particles initially close together in phase 
space separate exponentially. This continues until a time on the order of h~^, (33). After 
this time the phase of the particles may take on any values. We wish to ascertain the 
behavior of the particles for longer times. More precisely, if we consider an ensemble of 
particle trajectories, {Oj,Pj), then what are the moments, {pN — Po), {{pn — Pof')-, for large 
Nl (The brackets denote ensemble averages.) 

Firstly, let us consider the average force on the ensemble. 



^ {pN - po) 

N 



(44) 
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By shifting the origin of time by N (which leaves the difference equations invariant) we have 

/=^^^^, (45a) 
and by reflecting time using the symmetry, (23), we have 

/=^^^^. (45b) 

From (45a) and (45b) we see that / = 0. 

Turning to the second moment, we first define an acceleration. 

o-o=Pj+\-P3- (46) 
(In defining aj, we do not treat p as a periodic variable.) Then we have 

7V-1 N-l 
i=0 j=0 

We define a correlation function, C, by 

c(^^,i-J^={aiaj). (48) 

C has the following properties: 

C{e,k) = C{0,k) = Ck, (49a) 

which follows from the invariance of the system to a shift in the origin of time, and 

Ck = C_fc, (49b) 

since (aiaj) = {ajai). Because of (49a) the definition of may include an average over a 
single orbit, as well as an ensemble average. [We will usually use a subscript, N, to denote 
the time (i.e., the iteration number), while a subscript, k, will denote a time difference.] 
Rewriting (47) in terms of Ck and using properties, (49), we find 

Af-l 

{{pN-pof)=NCo+Y,'^iN-k)Ck. (50) 
fc=i 

For N ^ oo we have 

iiPN - Pof) = 2VN, (51) 
where the diffusion coefficient, V, is given by 

^ oo 

'D=^Co + J2Ck (52) 

k=l 

and we have assumed that Ck decays sufficiently rapidly so that the sum in (52) exists. 

In order to find Ck numerically we compute M orbits of length, L. Then Ck is approx- 
imated by 

L-k 

Cfe = (i-fc + l)-i ^(a^a^+fc), (53) 

e=o 

where the average is now over the M orbits. 
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Figure 8 shows examples of the correlation function for S = 0.23 and various values A. 
When A < As, Fig. 8(a), Ck is highly oscillatory, and decays slowly. As A is increased above 
the stochasticity threshold, the correlation time, kc, rapidly decreases. Fig. 8(b), becoming 
approximately for A > 1, Fig. 8(c). 

Note that in defining the ensemble of orbits we did not specify what po was, so that dif- 
ferent members of the ensemble could have different values of po- Because of the periodicity 
of the difference equations, (13), in p, we know that po (mod 2tt) is sufficient to determine 
pN — Po- Thus, the ensemble average provides an average of the initial conditions over an 
interval of 27r, which is of little consequence. Averaging along a given orbit (for A > Ag) 
has the same effect. Having performed this averaging, Ck is independent of p. 

The situation is much more complicated in the original system described by the Lorcntz 
force law, (2). In that case the problem is not periodic in r, since A has a slow dependence 
on r, (14). Now, we would be hiding this important dependence on r if we did not restrict 
the starting positions of the orbits, tq. Similarly, we could not perform the averaging along 
the orbit. Thus, Ck becomes a function of r, and it is more difficult to compute, since less 
averaging can be done. Finally, the diffusion coefficient would be more difficult to define 
because, although ((r^v — ?'o)^) might be proportional to N for some range of N, it certainly 
will not be proportional to N for large A'^, since the particles will be coming into regions 
where A and Ck are different. 

These problems arc all circumvented by using the difference equations. These allow a 
simple and rapid determination of Ck and V. The dependence of the diffusion coefficient 
on r will be recovered through its dependence on A (Sec. IX). 

VIII. THE DIFFUSION COEFFICIENT 

We now turn to the diffusion coefficient for the ions. We shall use (52) to define and 

calculate 2?. In the limit A ^ Ag wc can obtain an estimate for V based on the observation 
of the last section that the correlation time is short in this limit, viz., = for |fc| > [see 
Fig. 8(c)]. In this case we have T> = ^Co, with Cq = (oq), where the average is taken over 
an ensemble and down a trajectory. In the same spirit that wc calciilatcd the KS entropy, 
h, in the limit A Ag (Sec. VI), we shall assume that the ion trajectories are ergodic so 
that the ensemble and trajectory averages may be replaced by a phase space average. Thus, 
we obtain 

V = l{al) = l{[{v,-vo)-iu,-uo)r) 

= ^TT^A^{{cOSUi -|-COSWo)^) 
/■27r 

= 7r2^2(27r)-^ / cos^V# 
Jo 

= \^'A\ (54) 

When A is not large we expect the presence of islands to obstruct the diffusion of ions 
leading to a value of I? below (54). [This is reflected in the correlation function by the fact 
that for small fc, Ck>o < 0; sec Fig. 8(b).] Finally, for A < Ag diffusion is stopped. (Strictly 
speaking this only happens when 6 is not close to or ^, for then there are KAM surfaces 
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preventing the diffusion of ions. However, as disenssecl in Sec. V, the diffusion for A < Ag 
for S close to or 5 is extremely slow.) We therefore expect V to have the form, 

I?=i7rW(^), (55) 
where g{A) < 1, g{\A\ ^ 00) ^ f , and g{\A\ < As) = 0. 

We test these ideas out by numerically finding V, and hence g{A), for various values of 
6. We determine V by first finding Ck using (53) and then computing 

K 

I?= iCo + X^5feCfc, (56a) 
fe=i 

where qk is a windowing function given by 



5fe = <( 1 ^ f k-K' 

1 + cos TT- 



K + l-K' 



for < fc < K', 
for K' <k< K. 



(56b) 



The reason for introducing a windowing function into (56) is to suppress the effect of the 
oscillations seen in Fig. 8(a). Criteria for choosing K and K' are K — K' ^ 1 for A < Ag 
(in order to suppress the oscillations in Ck), K' ^ kc for A> Ag (so as to include the major 
contribution to V), and finally K L , where L is the length of the orbits used to calculate 
Ck, (53). This last condition arises because, from (53), the number of samples entering into 
Ck is M{L— k+1) which we wish to be large in order to suppress the statistical fluctuations 
in Cfe. 

Note that if we had computed D directly by measuring (pi — po)^ for M orbits, there 
would be only M samples contributions to V. Thus, the error in V would be 0{M^^^^). 
When computing V using Ck we effectively include all possible shorter orbits that make up 
the orbit po to p^. There are about L/kc independent orbits in an orbit of length, L. The 
error m V is then 0[{kjML)^/^] which can be much less than the error when computing 
V directly. 

Figure 9 shows g{A) for three value of S. Although g has the general form discussed 
earlier, there are some anomalies. Instead of approaching unity monotonically as A ^ 00, 
there is a tendency for g to oscillate about unity with a period of unity in A. This is 
particularly noticeable with S equal to 0.11 and 0.47. The oscillations arc about ±20% for 
1 < A < 2 and about ±10% for 20 < A < 22. Also there is a very strong peak in g{A) for 
6 = 0.47 and A = 0.55. Both these phenomena are due to the existence of "accelerator" 
modes, or island systems for which the ion, rather than returning to its original island, 
returns to an island displaced upwards or downwards in p by an integer multiple of tt. 

Some of these accelerator modes are studied in Appendix B. From (B4) with m = 1 
and n = we find that there exists a stable accelerating fixed point for S = 0.47 and 
0.53 < A < 0.5971. So for A = 0.55 there are orbits in the island system for which 
Pat ~ ±TrN. Although the orbits used to calculate V all lay outside this island system, they 
could wander close to the island and become temporarily trapped close to the island (for up 
to several hundred iterations). This results in a very slow decay of the correlation function 
(Fig. 10), and an anomalously high diffusion coefficient. 

As A is increased, a new set of accelerator modes appear at values of A equal to p±6, 
where p is an integer. The islands disappear by becoming unstable, shortly afterwards. This 
leads to the oscillatory behavior in V seen in Fig. 9. Similar oscillation have been observed 
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by Chirikov'' for the so-called standard mapping. He also associated these oscillations with 
the presence of accelerator modes. 

In our case there are probably no accelerator modes present in the cases for which g{A) 
is plotted in Fig. 9 (except the case previously discussed, i.e., S = 0.47, A = 0.55). The 
oscillations may however arise from the unstable fixed points that remain after the islands 
turn unstable. Such unstable fixed points may marginally increase the correlation time and 
hence the diffusion coeficient. As 5 is increased these fixed points become more unstable 
and the effect is diminished. Choosing A in the small windows where the accelerator modes 
exist would have given peaks in g{A) similar to the one at 5 = 0.47 and A = 0.55. We 
may argue that the peaks in g{A) that occur when the accelerator modes are present are 
unphysical. In Appendix B we establish that the window in A for which these modes exist 
is quite narrow. Now, as an ion is being accelerated or decelerated by one of these modes, 
p and r, (6a), will be changing. Hence, after a while we must recognize that A, which is a 
function of r, (14), will also change. This will "switch off" the accelerator modes after only 
a few iterations, and the correlation times will similarly be limited. The effect of the slow 
change of A (which was ignored in deriving the difference equations) will be to remove the 
peaks in g(A), leaving only the gentle decaying oscillations in g{A). 

For these reasons g{A) was fitted by the simplest smooth curve that reasonably fits the 
data. Figure 9 shows the functional form chosen, viz., 

5(A) = max (^—^,0j, (57) 

with As = J. This form of g will be checked in the next section by solving the resulting 
diffusion equation, and comparing the results with the solution of the exact equations of 
motion, (2). 



IX. DERIVATION OF A DIFFUSION EQUATION 

With A > As we expect the distribution, f{p, N), of a group of ions after N iterations 
to spread according to a diffusion equation. We will briefly outline the arguments that lead 
to this result. 

Suppose m is greater than the correlation time, kc, and let us look at the results of the 
mapping only at times, N, which are an integer multiple of m. Then the mapping, T, may 
be approximated by a Markovian process, which is described by a transition probability, 
P{Pi I P2, AT) dp2, which is the probability the a particle is found in the range, {p2, P2 + dp2), 
given that it was at pi, N iterations earlier. P must satisfy 

P{p,\p2,N)= j dpP{p,\p,k)P{p\p2,N-k) (58) 

for all k such that fc is a multiple of m and < k < N. 

For long times / will be spread over several periods in p. In that case we are not 
interested in the details of / that occur over the period, 2tt, in p. We retain the important 
physical effects when we average / and P over a period in p. (Precisely this averaging is 
carried out in computing Cfe and V in Sees. VII and VIII.) In that case P only depends on 
the difference, P2 — Pi, and we have 



P{p,\p2,N)=P{p2-pi,N). 



(59) 
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The evolution of the distribution of ions, /, is given by 

f{p,N) = j dpU{PuNr)P{p-p,,N-N,). (60) 

The moments of P, 

an{N) = J dpp"P{p,N), (61) 
are the same as the moments, ((p„ — po)"), studied in Sec. VII. Thus, we have 

ai{N)/N = 0, (62a) 
a2{N)/N = 2V, (62b) 

for N > m. 

Following the analysis of Chandrasekhar^^ and Wang and Uhlenbeck^^ we may use (58) 
to derive a Fokker-Planck equation for the evolution of P for N ^ kc- 

dP 5^ 

Substituting for / from (60) gives 

df d d 

where we have used the result that V is independent of p in order to write the diffusion 

equation in its more usual form. Incidentally, although wc derived (64) for A > Ag it clearly 
gives satisfactory results for A < Ag, because g and T> are then zero, and / does not evolve. 

We now undo the normalizations made earlier. Writing TV = i/27r and converting p to 
r using (6a) we have 

dt rd/2TT\dr) dr^' ^ ' 

where dp/dr = {r^ — z/^)^/^/r = Hi^''' (r) / Hi^\r) for r > + (^i/)^^^. In writing (65) in 

this form we have replaced the first derivative operator in (64) by the cylindrical divergence 
operator, in order to ensure conservation of particles. (The normalization of / is such that 
J J dr 2nrf = 1. Wc will suppress the dependence of / on the velocity parallel to B, . 
until we consider the effects of collisions in Sees. X-XII.) Note that T> is now a function of 
r through A, i.e., V = V[A{r)], so that we must justify our commuting T> with the second 
derivative operator between (63) and (64). The argument for having V where it is in (65) is 
that in the steady state / should be a constant in the stochastic region. This is a result of 
the approximate ergodicity of the motion, when described by the exact equations of motion, 
(2), and was observed in the simulation described in I. Wc shall presently compare the 
results obtained by the diffusion equation with that simulation. Rewriting (65) using (55) 
and (14) we obtain 

dt r dr^ dr ^^^^ 



where 
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for r>u. We extend the definition of Z) to r < as follows: 



{D{v), for V — y/a <r<u, 

D{v) [r~{v- 2yfa)f/a, for v - <r <u-^/a, 
0, iov r <y - 2\/a. 

When numerically computing D we use the approximate forms for H and if', 



(67b) 



1/2 ( 



-1/4 



for r > 1/ + {\vY^^ and 



if, 



(68a) 



(68b) 



(68c) 



for 1/ < r < di')^/^ and similarly for H' . We will show the typical form of D in Sec. XL 

Because of the restrictions on r in the derivation of the difference equations (Sec. II), 
we must justify (67) for r < u + (i;/)^/^. In this range, the choice of the form of D is 
motivated by the following considerations: Firstly, in the region, v — y/a <r <v + {\vY^^ 
trapping^ plays a dominant role. This trapping is similar to the trapping of a particle by a 
wave in the absence of a magnetic field and its effect is to cause a rapid mixing of / in this 
region. This is modeled by (67) since D is approximately constant and has its maximum 
value in this region. Secondly, we find that particles are fed quite slowly into the trapping 
region from below and this is given by the form of D iov u — 2y/a < r < u — y/a. Finally, 
the motion is coherent for r < v — 2y/a in which case D is 0. 

In order to check (66) and (67) we compare the results obtained using these equations 
with the results of the simulation described in I. In the latter case the orbits of 50 particles 
were found by integrating the exact equations of motion, (2). Because of the small number of 
particles involved the comparison is most easily made if we solve (66) using a simple Monte 
Carlo method described in Appendix C. It might seem as though wc are back-tracking here, 
by going back to difference equations. However, the Monte Carlo method for (66) differs from 
the original equations, (13), in two respects: firstly, the randomness is explicitly inserted; 
and, secondly, we allow for the variation of D with r. This latter aspect of the solution we 
present here allows us to justify the placement of D with respect to the derivative operators 
in (66). Figures 11 and 12 show the comparison. In each case we follow the motion of 
TV = 50 particles with initial perpendicular velocities, ro = 23, and with a = 20, ly = 30.23. 
The general features of these figures are discussed in I. For our purposes we see that there 
is good agreement between the two simulations. In particular, note the agreement in the 
short time behavior, Fig. 11(c) and (d), and in the averaged distribution function, Fig. 12. 
These two points respectively confirm the form taken for D for r < in (67b) and the 
placement of D with respect to the derivative operators in (66). We conclude that (66) with 
D given by (67) accurately models the diffusion of ions in the presence of a perpendicularly 
propagating electrostatic wave. 

We may write D in dimensional form by multiplying (67) by to give 
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4' 



klvl 

D{v±) = I D{vph), for Wph - vtr <v±< Uph, (69) 

■D(Wph)b_L - {Vph - 2vtr)]'^/v^^, for Wph - 2vtr <V±< Vph - Vtr, 

.0, for vi_ < Vph - 2vtr, 

where g{A), v, and r axe given by (57), (14), and (3), Vp^ = ui/kj_, and Vt^ is the trapping 
velocity, 

Vtr = {q^Eo/m,k^y/\ (70) 

In (69) we have included an additional factor, 7, which is defined as the fraction of time 
that the ion orbits spend in the region of lower hybrid waves. This reflects the fact that the 
ion can only diffuse a fraction, 7, of the time. In order to include the spatial variation of the 
lower hybrid waves in such a simple way, two conditions must be satisfied. The ions must 
spend many cyclotron periods in the lower hybrid waves, so that the difference equations are 
applicable, and the motion is described by (64). The opposite limit, where the ion spends 
a small fraction of a cyclotron period in the wave, has been considered by Lazarro^^ in the 
context of heating a fast beam of ions. Secondly, we must be able to model the lower hybrid 
wave as a product of a single k component and a square-wave envelope. The first condition 
is usually satisfied in cases of practical interest and, because the lower hybrid waves travel 
in well-defined rays,^^ the second condition is also normally satisfied. In cases where the 
envelope of the lower hybrid waves is not a square wave, the definition of D, (69), may be 
replaced by the average of (69), without the factor, 7, over the ion's orbit. For a circulating 
ion in a tokamak which covers a magnetic surface ergodically, 7 is approximately the ratio 
of the area of the intersection of the lower hybrid ray with the magnetic surface to the total 
area of the magnetic surface. This ignores the finite extent of the lower hybrid wave in 
the perpendicular direction. When this is included, 7 is reduced since the cyclotron orbit 
must be completely in the lower hybrid ray for the difference equations, (13), to hold. The 
diffusion equation for the ions written in unnormalized form becomes 

dt 



v±_ ov± ov± 



(71) 



In the next sections we study the properties of this equation when we add a Fokker-Planck 
collision term. 



X. INCLUSION OF COLLISIONS 

We have seen that the ions can irreversibly exchange energy with the wave. However, 
because u}/k± is usually several times the ion thermal speed, only the tail ions are affected. 
Bulk heating only takes place when these ions collide with the background. We investigate 
this process by including a Fokker-Planck collision term in (71) to give 

df 15 „ 5 , /df\ , , 

-vi.D — ]+[^\ . (72a) 



dt v± dv± dv± \ dt 

Since only a relatively few tail particles are affected by the wave we can neglect the tail-tail 
collisions when determining {df/dt)c- We can thus linearize this term by assuming the the 
background distributions of ions and electrons are Maxwellians. We then obtain^^ 

{df/dt)e = -yZyr/^, (72b) 
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where /? is a species label {i or e), 

+ rufj 4 

and expressions for the collision frequencies, v, can be found in Ref. 18. 

Wc next simplify (72) by reducing it to a problem in only one velocity dimension, 
by integrating (72a) over I'll . The procedure we follow for doing this is very similar to that 
employed by Fisch^^'^° to derive a one- dimensional Fokker-Planck operator in W|| in order 

to study rf-driven currents. When performing the integral over t;|| the term involving J^^^ 

drops out. We then need only to evaluate — / dv\\ . In order to do this we make three 
assumptions: Firstly, we assume that / {— fi) has the form, 

f{v^,v\\) = F{v±){2TrTi/mi)-'^^eM-h"'ivl/Ti); (73) 

i.e., the dependence of / on W|| is that of a MaxwcUian. Secondly, wc assume that v'^ ^ 
Ti/rrii. This allows us to replace v by v± in the expressions for the collisions frequencies, so 
that they drop out of the integral. Lastly, we assume that the background temperatures are 
all equal so that = Ti. (Without this assumption additional sources and sinks of energy 
are required to maintain an equilibrium.) Then we have 

- / dv\\ fl^ = Cn dF/dv^ + {miV^/Ti)C0F + 0[Ti/{mivl)], (74) 



where 



and we have used^^ 



C0 = \v\'^vl + \v'l^TJm„ (75) 



TOq + TO/3 * Tfj 2 II 



Substituting these results into the integral of (72a) over vn gives 



dF _ \ d 

dt v±_ dv± 



OV± \OV± I' 



Li 



(77) 



where C = ^^'^ ^ti = Ti/rrii. If D = 0, we recover a Maxwellian for F. Thus, we 

shall use (77) to solve for F for all v± even though it was derived only in the high velocity 
limit. We expect that this will not entail much additional error as long as D is finite only 
where {v±/vti)^ is large. 
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XI. NUMERICAL SOLUTION OF FOKKER-PLANCK EQUATION 

A number of approximations need to bo made in the derivation of tlie one-dimension- 
al Fokker-Planck equation, (77). These leave the accuracy of (77) open to question. In 
particular the factorization of / by (73) is not easy to rigourously justify. In this section we 
compare the numerical solutions of the one-dimensional and two-dimensional Fokker-Planck 
equations, (77) and (72). A similar, but more extensive, comparison has been carried out^^ 
for the parallel one-dimensional Fokker-Planck equation for the electrons with a parallel 
quasilinear diffusion term, derived by Fisch.^^'^° For reasons that we will discuss wc; cixpcct 
better agreement in the present case. Therefore, we content ourselves with presenting here 
the solution for one particular case. 

Consider a hydrogen plasma with ion density, no = 10^° m~^, magnetic field, Bq = 
4 T, and temperature, Tg = = 2 keV. Let the lower hybrid wave be described by 
oj/{k±Vti) = 3.5, uj/ojih = 1.3, and Eq = 10® V/m. The frequency of such a wave is 
2.13 GHz. Its parallel wavenumber may be found by solving the dispersion relation -|- 
K±k'j^ — 3(ti'pju|j/u;^)fcj^ = 0, where and K± arc elements of the cold plasma dielectric 
tensor. This gives u)/{k^\Vte) = 5.25 and a parallel index of refraction, ny = 3.06. This 
wave represents a typical Fourier component of a lower hybrid wave near the center of the 
tokamak, but before it has reached the point of thermal wave conversion. Note that, because 
of the high value of ui/{k\^Vte), we expect electron Landau damping to be negligible. The 
normalized parameters are then found from (3) tohe v = 35, a = 5.71. We take the extent 
of the lower hybrid ray in the parallel direction to be Az = 0.4 m, the major radius of the 
tokamak to be i? = 2 m, and the angular extent of the wave in the poloidal direction to be 
= 7r/4. The geometrical factor, 7, appearing in (69) is given by 

^ Az Ae _ 1 
^"^^ 2^2^ ~ 251' ^ ' 

[As discussed in Sec. IX, this value of 7 should be slightly reduced to account for the finite 
perpendicular extent of the waves. However, (78) is sufficiently accurate for illustrative 
purposes.] Substituting these values into D, (69), at v± = ui/k± gives 

£>(a;/fci) = 1.5 4t/f , (79) 

3 of a collis 

ion collision frequency for thermal particles 



where we have written D in terms of a collisional diffusion coefficient v^^Uq^ and is the 



= (80) 

We take the Coulomb logarithm to be Xu = Aje = 15. With our parameters = 4.31 x 
10^ s~^. The full shape of D is given in Fig. 13. This is used in the numerical integrations 
of the Fokker-Planck equations. 

We now solve the one-dimensional and two-dimensional Fokker-Planck equations, (77) 
and (72), numerically. Att = both / and F are taken to be Maxwellians with temperature, 
Tj. The diffusion coefficient is turned on with a linear ramp function between t = and 
t = {i^Q^)~^ and is kept constant thereafter. 

Figure 14 shows the steady-state distributions. [In plotting the perpendicular distribu- 
tion function in the two-dimensional case we have defined 

F{v^) = J dv\if{v^,v^\). (81) 
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Note that this is consistent with (73).] Although the dependence of / on ?;|| is not that 
of a MaxwelUan, the two results for F{v±_) are in close agreement. The major difference 
is that the tail of the distribution function is characterized by a temperature of Tj in the 
one-dimensional case, but by about 1.16 Tj in the two-dimensional case. 
The power dissipated by the wave is defined by 

Pd = -2-nnomi j v'iDdF/dv±dv±, (82) 

while the power lost by the test distribution to the background distributions is 

= norui j V j^/'^ dv (83) 

in the two-dimensional case and 

Pe = 27rnomi j vlC[dF/dv^ + {v^/vl)F] dv^ (84) 

in the one-dimensional case. These are plotted as a function of time in Fig. 15. The most 
important quantity we wish to know is the steady-state power dissipation, P. As t — > oo, 
we have P = P^ = Pc- We see that there is close agreement between the one-dimensional 
and two-dimensional values for P. (They differ by about 3%.) However, the solution to 
the one-dimensional equation reaches a steady state about 2 times faster than that of the 
two-dimensional equation. From Fig. 15(a) we see that Pd nearly reaches its steady-state 
value at i = Td « ^/^o^ ^ 1.4 ms, whereas Pc (and / also) reaches a steady state on a time 
scale, Tc, which is about or 14 ms. 

The agreement between the one- and two-dimensional treatments in the power dissi- 
pated is much better than observed in the case of the parallel one-dimensional Fokker-Planck 
equation for the electrons. This is so because the pitch-angle scattering terms, which are 
the principle reason that / deviates from the form assumed in (73) are approximately 4 
times less important in the case considered here. Two effects contribute to this factor of 
4: In our case electrons do not cause appreciable pitch-angle scattering; whereas, in the 
other case, ions and electrons contribute equally. The other effect is geometrical in origin 
and arises because the perpendicular plane has two degrees of freedom, as opposed to the 
parallel direction, which has only one. 

We conclude that the one-dimensional Fokker-Planck equation, (77), gives accurate 

results for the steady state distribution and for the power dissipated in the steady state, 
P. The agreement on the time scales is worse, although we expect the one-dimensional 
equations to give qualitatively correct results here. 



XII. ANALYTIC SOLUTION OF THE ONE-DIMENSIONAL EQUATION 

We now discuss the analytical solution of the one-dimensional Fokker Planck equation, 
(77). We shall be mostly concerned with the steady-state solution, since this is of the most 
practical interest, and accurately approximates that givc^n by (72). We will bric^fly consider 
the time scales given by (77). The methods we use and results we obtain are very similar to 
those of Fisch;^^'^" who treats the parallel one-dimensional Fokker-Planck equation. Most 
of the analysis will be carried out for general collision term, C. We will however sometimes 
need to know the functional dependence of C on v±. We shall therefore give the form of 



-22- 



C in the limit vu <C <C Vte- Using the forms for v^^^^ and v^^^ given in Ref. 18 and 
assuming that =Ti, which is required in the derivation of (77), and Aje = we have 

where is given by (80) and 



9Zi f rrii^ ^^"^ 



or 

vo_ 

Vti 



(86a) 



'6.23, for hydrogen, 

6.99, for deuterium, (86b) 
^7.48, for tritium. 

The contributions of Ci and Ce to C are given by the first and second terms in parentheses 
in (85). 

The steady-state solution to (77) is^° 

F(t ^ 00) = Fo exp - YVdIc '^''^ ' ^^^""^ 
where Fq is a normalization constant determined by 

j 2'!TV±Fdv± = 1. (87b) 
If D is only non-zero for v'^ ^ v^^ then 



Fo « (27r4)-^ (88) 
The power dissipated in the steady state is given by (82), 

P=^-^ fr^^Fdv^. (89) 



1 + D/C 

This power dissipation is divided between the electron and ion backgrounds according to 

2Trnomi f C^Dv^^ 

where is the power deposited in the background species, 0. 

Equation (89) may be evaluated in two limits. Firstly, in the collision dominated regime, 
for which D -C C for all v±, F is approximately a Maxwellian and (89) becomes 

1 
2«? 



noTUi 



P=^ I Dvi exp ( --^ ) dv±. (91) 



If D{v±) = for v± < vi and D{v±) ~ Dq (a constant) for vi < v± < vi + v^^/vi, the 
integral may be approximated by Laplace's method to give 

P « nom^ivi/vufDo exp{-^vl/vl). (92) 

If vf -C fo, nearly all of this goes to ion heating. 
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The second limit we shall consider is the "plateau" limit, in which a plateau forms on 
F. Suppose -D = for < v\ or i>j_ > vi (with vi > v\) and D ^ C otherwise. If D also 
satisfies 

dw± « 1, (93) 

then F is approximately constant between vi and V2- Assuming Vi » v^^, so that Fq is 
given by (88), (89) becomes 

Substituting for C from (85) we have 



nomiiyl/\)u[{v2 - wi) + - 'f]')/wo] cxp(-it'i/tf,). (95) 



where the first and second terms in brackets are the contributions of P to the ion and 
electron heating respectively. 

The time scales may be readily computed for this case. Since the time scales derived 
from the one-dimensional equation are not very accurate, we shall merely quote the results. 
The derivation closely follows that given in Ref. 19. The time for the distribution to flatten 
between vi and V2 due to D is 

Td = {v2-vi)yD. (96) 
If Vi{v2 — Vi) » Wjj, the time scale for particles to collisionally flll in the plateau is 

re = U^l-vf)vy[vlC{v,)], (97) 
or substituting for C, (85), with Vi <^ Vq 

rc = U4-v!)v^/{v!A^'). (98) 

The time scale for the saturation of Pc is Tc, while is the time at which Pa nearly reaches 
its steady-state value. For t > Td the power required for the wave to accelerate the ions 
into the plateau decreases, while the power to maintain them there increases. Thus, Pa is 

roughly constant for t > Td- 

Many easels of interest will not be in the collisional or plateau limits. However, we may 
still obtain approximate results by judiciously taking one or other limit. For instance, if 
we consider the case treated in the previous section, we see that a plateau has formed over 
part of the range where D is finite. Fig. 14(b). The plateau limit has some validity in this 
case since D cuts off quite sharply at the low velocity end. Fig. 13, and so the height of the 
plateau is fairly well determined. We take v\ to be given by 

D{vi)=C{vi) (99a) 

and V2 by 



V2 



^dv^ = l. (99b) 



(Thus, at -F is a factor, e, lower than the plateau height.) These give V\ = 2>.lvti and 
V2 = 5.9 Vti, approximately. We then obtain 

P « 5 X 10^2 mmivlivll". (100) 
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About 70% of this goes to heating the ions. Equation (100) is fairly close to the observed 
value for which the numerial factor is 3.8 x 10~^, Fig. 15(a). By far the most sensitive factor 
in determing P is since the height of the plateau depends exponentially on vi. More 

accurate estimates of the plateau height are possible. However, if a more accurate result 
that (100) is required, it is probably easiest just to evaluate (87) and (89) numerically. 

With these values of vi and V2 the time scale become Td = S/i/^* [taking D to be given 
by (79)] and Tc = l^jv^^ . These agree quite well with the results of the numerical solution 
of (77) given in Fig. 15(b). However, as we noted these times (and, in particular, Tc) differ 
from the times given by the two-dimensional Fokker-Planck equation. 

XIII. DISCUSSION 

We have computed the rate of diffusion of ions in perpendicular velocity space in the 
presence of a uniform magnetic field and a perpendicularly propagating electrostatic wave 
for which a; ^ 57^. The diffusion coefficient is given by (69). From the results of 1, it is seen 
that this diffusion coefficient also applies to diffusion in a lower hybrid wave (for which k\\ is 
finite, but fc|| <C fcx) in a weakly inhomogeneous magnetic field. The effect of nonuniformity 
of the lower hybrid wave is modeled by introducing the geometrical factor, 7, into (69). 

In order to determine the heating rate, collisions are included, and a one-dimensional 
Fokker-Planck equation, (77), is derived. Prom this the steady-state power dissipation, P, 
may be found, (89). Provided lo / {k±^vti) becomes sufficiently small (on the order of 4) this 
may result in strong damping of the lower hybrid wave and efficient heating of the ions. 
(See the example in Sec. XI.) 

These results may be used in a simiilation of lower hybrid heating with a transport 
code, provided the time scale for the evolution of the plasma exceeds Tg. The quantities, Pg 
and Pi, (90), which represent the powers deposited per unit volume, would then appear as 
source terms in the equations for the electrons and ions. The wave would obey an evolution 
equation 

Vgr dW{r, t)/dr = -P(r, i)/7, (101) 

where W is the energy density of the wave, r is the minor radius of the tokamak, and Vgr 
is the radial group velocity. In this equation t is a parameter. P and W would vary with 
time as the background distributions evolve. 

Lastly we consider the question of particle transport due to stochastic heating. Firstly, 
there is the change in the transport coefficients due to the altered distribution function of 
the ions. Fig. 14. Since the heating initially takes place in the perpendicular direction, this 
would lead to an increase in the number of trapped ions, which may have a deleterious 
effect on the transport. Secondly, the equations of motion, (2), give a shift in the guiding 
centers of the ions. The amount of this shift is most easily calculated from the Hamiltonian, 
(Al), since —I2 is the x guiding center position. (There is no shift in the y guiding center.) 
Thus, we have — A72 = AIi/v (after averaging over a cyclotron period) or, undoing the 
normalizations 

Ax = (fc_L/w) A/xk X B, (102) 

where fj, is the magnetic moment, ^miV^/{qiBo). When considering heating by the injected 
lower hybrid wave (rather than parametric decay products) k is mainly in the radial direction 
and so Ax is mostly in the poloidal direction. So this effect does not lead to any outward 
transport of ions. 
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APPENDIX A: DIFFERENCE EQUATIONS FROM THE HAMILTONIAN 

In this appendix we derive the difference equations, (13), from the Hamiltonian for the 
ion. Using this approach it is easy to establish what approximations are made. In I it was 
shown that the motion of an ion in the fields, (1), is governed by the Hamiltonian, 

H = Ii + uh — asin[(2/i)^/^ sin — 1^2] 

00 

= h + i^l2-a J2 ^m[(2/i)'/']sin(mu;i-u;2), (Al) 

m=— 00 

where the normalizations arc the same as used in (2), Jm is the Bessel function of order m, 
y = (2/1)^/^ sin wi, and x = —I2 — (2/i)^/^ cos Wi. The action variable, Ji, is related to the 
Larmor radius, r (3c), by Ii = ^r^. We now perform a series of transformations, similar to 
those made in Appendix C of I, to bring the sum in (Al) into the form of an infinite sum 
of the product of two trigonometric terms. Hamilton's equations of motion may then be 
written in terms of delta functions. 

We first transform (Al) to hatted variables using the generating function, 

F2 = Ii{nwi - W2) + I2W2, (A2) 

where n is given by (15). Equation (Al) then becomes 

H = -5ii + 1//2 - a ^ J„[(2n/i)i/2] sin 

m 

Since Wi is slowly varying, whereas W2 is still rapidly varying (assuming n ^ 1) the dominant 
contributions to the sum in (A3) come from m ~ n. Recognizing this we replace m by n + A; 
in the sum to obtain 

N 

H = -5h+ul2-a J2 ^+fe[(2n7i)V2]sin 

k=-N 

where N is some integer satisfying 1 N n. [It is convenient to choose N w (i;^)^/^.] 
As a result of this inequality we approximate the factor (1 + k/n) by unity. We next use 
the large argument expansion of the Bessel function 

J„(r) = (2/7r)V2(,2_ ^2^-1/4 

X cos[(r^ - m^y/^ - m. cos^^{m/r) - 7r/4] (A5) 

[Ref. 3, Eq. (9.3.3)], which is valid for r > m + (im)^/"^. In order to put the Hamiltonian 
into the required form we need to be able to do three things: approximate all the Bessel 



m 



-Wi 



— I W2 

n, J 



(A3) 



-W2 



(A4) 
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functions in (A4) using (A5), which requires r > (n + A'') + (J^nY^^; replace the amplitude 
factor (r^ — to^)~^/^ by (r^ — i^^)^^/^, which requires r"^ ~ -n? ^ Nn; and to Taylor-expand 
the argument of the cosine in (A5) in k about k = S, keeping only the first two terms, which 
requires r"^ — ^ N^. With N « (^i/)^/^ all these inequalities become r — u ^ (ij/)^/^. 
The Bessel functions in (A4) then become 

Jn+kir) = (2/7r)i/2(r2 - ;.2)-i/4 cos(i? - k^ + S<j)), (A6) 

where r = (2n/i)^/^, R = g{r) = (r^ — j/^)^/^ — cos~^(z//r) — 7r/4, and (j) = cos~^{v/r). 
We shall regard the r's appearing in the amplitude factor and in the definition of (j) as 
being parameters. This requires that neither of these quantities change significantly when 
r changes by a period of the Bessel function, 27r{dR/dr)~^ . The conditions on r and v 
already assumed are sufficient to guarantee this. In addition, if a exceeds the stochasticity 
threshold, r may change by more than the period of the Bessel function, in which case we 
must re-examine our assumptions about the constancy of the amplitude factor and (j). We 
shall return to this point after we have derived the difference equations. 

We now transform to a new set of variables, distinguished by tildes, using the generating 
function, 

F3 = -G{ii)wi - i2W2, (A7) 

where G(/i) = g[{2niiY^^]. This gives 

/i = G(/i), wi=G'{h)wi. (A8) 

We solve both of these locally by expanding Ii about some point Iio and keeping only the 
terms involving G'{I\q). Thus, we obtain 

h = h/G'{iw) + const., wi = G'{iw)wi. (A9) 

This is consistent with the approximations already made. Writing G'(/io) in terms of the 
Larmor radius, r (which, like /lo, we regard as constant), wc have 

G'(Jio) = g'ir) dr/dii « (r^ - u^^^u/r'^ = Q. (AlO) 

The Hamiltonian becomes 

H = --ii + iyi2 
Us 

Finally, we perform the scaling transformation, 

M = QH, (A12a) 

Ji =/i, - Qwi, (A12b) 

J2 = Qvh, i)2 = W2IV, (A12c) 
so that the Hamiltonian is given by 

00 

M = -8Jx^ J2 - A ^ cos(Ji -A:(?;' + (50)sin(V'i +A:V'2), (A13) 

fc= — 00 
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whcrc A is given by (11) and wo have approximated vjn by unity to give the last term in 
the argument to the sine. We have extended the hmits of the sum to infinity; the additional 
terms introduced are non-resonant, and so do not have much effect. 

Hamiltons's equations now give ^"2 = 1 or ^/)2 = t and 

V>i = -5 + >iy^ sin(Ji - A:(/>+ (5(/))sin(V'i + A:t), (Al4a) 

j\= cos( Ji -k4>-\- 5(f>) cosfV)! + kt). (A14b) 

Defining new variables, 



we find 



^ = mr — tpi, p = Ji + i^TT, (Al5a) 
u = J — p, V = 'y + p, (A15b) 

u = 5 - A cos[w - (tt - 0)5 - k{t + (f>)] 

oo 

= S-2TrAcos[v-{Tr-(j))d] ^ 6{t + cj) - 2Trj), (A16a) 

j=-oo 

i) = 5 + Ay^ cos[-u - (tt - (f))5 + kit - (f>)] 

oo 

= 6 + 2TrAcos[u+{n-(j))6] ^ 5{t - cj) - 2Trj). (A16b) 

j=-oo 

We have used the notation, d, to denote the Dirac delta function (to distinguish it from the 
variable, 6) and we have used the identity, 

oo oo 

coskt = 2TT 5{t-2-Kj). (A17) 

fc— — oo J — — oo 

Defining Uj = u{t = 2iTj — tt) and similarly for v, we may obtain 

Uj+i — Uj = 2^5 — 27rAcoswj, (A18a) 

Vj+i — Vj = 2^6 + 2wAcosuj+i. (A18b) 

It is readily established that p as given by (A15a) agrees with (6a), and that ^ a,txp2 = 27rj— tt 
is equal to 1102 at wi = 2iTj — tt and so is equal to 9j as defined in (6b). Thus, the difference 
equations given in (AI8) agree with those derived in Sec. II, (13). The final condition on 
the validity of (AI8) is that the change in p in one iteration is insufficient to cause an 
appreciable change in A through r. (This then justifies taking A to be a constant.) The 
maximum change in p in one iteration is on the order of A. Therefore the fractional change 
in A due to this change in p is 

d\ogA_dA dAdr av 
dp ^ 'd^ ^ ^d^^ {r^ - 1/2)5/4 ' ^ > 
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whcrc, in differentiating A we have only considered its dependence on r through (r^ — v^Y/^. 
Demanding that (A19) be small gives 

a « - z.2)5/2/^ A « - J.2)3/2/^2_ ^^20) 

The other assumptions made are 1 and r — 

APPENDIX B: FIRST-ORDER ACCELERATOR MODES 

We consider here the simplest of the accelerator modes. Such modes are island systems 
around first-order fixed points in the (u, v) plane; i.e., points for which ui — uo = —27m and 
vi — vo = 27rm, where m and n are integers. The amount by which p increases per iteration 
is STT, where s = m + n. From (13) the fixed points are given by 

6 + Acosuo = m, S — Acosvq = —n. (Bl) 

In order for an island to exist, the fixed point, (uq^vo), must be elliptic. The stability of 
{uo,vo) is determined by the eigenvalues of the linearized mapping. From (35) we have 

Jo=(-2f/ 1-4^)' (^2) 

where U = it A sin uq and V = 7rAsint;o. The eigenvalues are given by 

A = 1 - 2UV ± 2(C/V^ - UVy/'^. (B3) 

The fixed point is elliptic if A is complex, which requires < UV < 1. Substituting for 
sin Wo and sinuo gives 

A^-{m- Sf > 0, (B4a) 

^2 - (n + (5)2 > 0, (B4b) 

[A^ - (m - Sf] [A^ -{n + Sf] < l/7r^ (B4c) 

as the conditions for first-order accelerator modes. When \A\ just exceeds max(|m — 
|n-|-(5|) then (B4) is satisfied. At this point elliptic and hyperbolic fixed points appear 
(with different signs for cosuq or cosfo). As |^| is increased further (B4c) fails and the 
elliptic fixed point turns into a hyperbolic fixed point with reflection. Thus, (B4) defines a 
range, A( < \A\ < A„, which is a function of 6, m, and n, in which first-order elliptic fixed 
points exist. 

Figure 16 shows the behavior of the mode with m — 1 and n = for S = ^. From (B4) 
wc sec that for these parameters Ai = ^ and Au = {j + 7r^2)i/2 — 0.5927. At A ~ 0.55, 
Fig. 16(a), the first-order island is visible. This is surrounded by higher-order (4th and 
32nd) accelerator modes. When A just exceeds A^, Fig. 16(b), the center of the first-order 
island becomes unstable, becoming a hyperbolic fixed point, separating two second-order 
islands. However, although this point is unstable, it is confined to the island system, since 
there are still eigencurves with the original topology surrounding the system of two islands 
and one hyperbolic fixed point. Thus, a point started close to the hyperbolic point will 
be forever accelerated. These eigencurves are soon destroyed, however, allowing the escape 
of particles starting near the fixed point. In Fig. 16(c) we see that the particle spends on 
the order of 1000 iterations encircling the second-order islands before escaping. As A is 
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incroasccl further, the second-order islands rapidly decrease in size, becoming unstable when 
A is slightly greater that 0.601. 

From Fig. 16 we see that we should study accelerator modes of many different orders in 
order to have a complete understanding of their effect on the diffusion of particles. However, 
while it is possible to thoroughly study first-order fixed points, just cataloguing the higher- 
order fixed points becomes a monumental task for A> 1. For example, there are over 100 
second-order fixed points present for ^4 = 1. We will therefore pursue only the first-order 
fixed points. 

In Fig. 17 we show the range in |A| defined by (B4) as a function of 6 with s = 0, 1, 
and 2. (For each value of s, all values of m and n satisfying m + n = s are considered.) The 
pattern evident in these figures is repeated for higher values of s. From (B4) and Fig. 17 
we see that an accelerator mode, with a given value of s, first appears at 

f|s/2| -t- 1^1 , for seven, 

Ae=\ (B5) 
\|s/2| + i-|5|, for s odd. ^ ' 

(Recall that 6 is defined such that \S\ < i.) 

The accelerator mode disappears at a slightly higher amplitude, but periodically reap- 
pears (with different m and n, but the same s) at 

{|s/2| ± \5\ +p, for s even, 

(B6) 
|s/2|±(i-|^|)+p, for s odd, ^ ' 

where p is a positive integer. For a given s, the gap in which the accelerator mode exists, 
AA = Au — Ae, is widest when A^ = \s/2\, 5 = for s even and \5\ = 5 for s odd. In this 
case we have 

A„=(-H-^) . (B7) 

For s = (i.e., no acceleration), we have A^ = Q and A^ = I/tt. This is the range in which 
the central fixed points of the islands in Fig. 3(a) are stable. If \s\ ^ 1 then 

When \A\ greatly exceeds the minimum threshold for a given s, so the A^ is given by (B6) 
with 1, we may find an approximate solution for A A from (B4): 

for s = 0, 



, . ,0 , for s 7^ 0. 
, An^sAf 

Since the number of first-order accelerator modes increases with increasing A (sec Fig. 
17), we are lead to ask how the fraction of phase space occupied by accelerator modes varies 
with increasing A. We try to find an upper bound on this quantity for the first-order modes. 
We estimate the size of an accelerator mode by noting that the maximum radius of the island 
is given by the distance between the elliptic and hyperbolic fixed points. We consider here 
the pair of fixed points that are born together and we assume that the islands are roughly 
circular. From (Bl) this distance is about [{A — Ai) / Ai]^/"^ (ignoring numerical factors). 
This is maximaized by letting A = A^- The maximum area of the island is then AA/A^. 
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For a given A, modes with s from 1 to about 2A can be present; thus, the total area of the 
accelerator modes is on the order of 

E^~Ei.~^=r^ (BIO) 

s— 1 s— 1 

Here we have used the form of for p 1 and s 0, (B9). If (5 = or i then there is 
an additional contribution of order from the mode that just appeared for the first time 
by satisfying (B5); in this case AA is given by (B8). This then dominates over (BIO). In 
either case the area of the accelerator modes decreases with A. 

However, the modes tend to line up close to cosw = and cosu = 0, since (B4a), 
say, may be well satisfied, when (B4b) becomes true;. We should investigate whether the 
modes can form a barrier inhibiting the diffusion of particles. The sum of the widths of the 
accelerator modes is 

2A 
s=l 

which also decreases with A. 

We conclude that the effect of the first-order accelerator modes decreases with increasing 
A. This leaves open the possibility that the effect of all the accelerator modes may not 
decrease. We have, however, observed no evidence of this. 



APPENDIX C: MONTE CARLO SOLUTION 

We describe here the solution of (66) by a Monte Carlo method. We first define a new 
velocity coordinate, ji = ^r^, so that (66) becomes 

(The variable fi is proportional to the magnetic moment and is the same as the action, I\, 
used in Appendix A.) We solve (CI) by dividing each cyclotron period up into M equal time 
intervals of 27r/M. The distribution function is recorded by the positions of N particles, 
/ij^j, where the subscript, i, refers to the particle number and j to the time step, j = Mt/2-K. 
The distribution function, /, is found by converting fi back to r and letting 

where SN is the number of particles with velocity between r — 6r/2 and r + 6r/2. We 
advance / by the difference equation, 

{+S~^, with probability |, 
—S , with probability 2- 

[Do not confuse 6^ with the parameter, 6, defined in (15).] In the limit, M — > oo, both 
and 6~ are given by 5{ii), where 

5(/i) = {27Tr^D/My^\ (C4) 

However, there are corrections to the expressions of and of order which must be 
included if we are to correctly model (66). [If we do not include these corrections we end up 



-31- 



solving the equation, df/dt — {d^ / d ii^)r'^ D f ] Wc derive these corrections by demanding 
that the steady-state solution given by the (C3) agrees with the steady-state solution to 
(CI), viz., / = const. Letting / be a constant (which we take to be unity), the flux of 
particles, S, at is 

S{^i) = \[5+{^i-e+)-5-{^i + e-)], (C5) 

where 

e+ =5+(/i-e+), e- = 5-{iJi + e-). (C6) 

Equation (C5) just expresses the fact that half the particles between — e+ and fi will pass 

the point ^ is one time step, and similarly for the particles between and ii + e~ . Setting 
S{^) to zero gives 

5-{^l + e) = 5+{^i-e) = e. (C7) 

Putting e = (5(/x) we recover the correct result for (5+ and 5~ in the limit, M — > oo. This 
gives 

5-[^i + 5{p)]^5+[^,-5{p)] = 5{^i). (C8) 

Equation (C8) has a very simple graphical representation which is shown in Fig. 18. In the 
limits, M — > oo and J ^ 0, we can make a Taylor series expansion about to give 

5+{ii)-5-{iJL) = d{5'^)/di^. (C9) 

This acts as a friction force which appears in the diffusion equation, (CI), when we write it 
in the form. 



dt dfj, 



(CIO) 



Prom Fig. 18 we see that 5+ or 6~ becomes double-valued if \dS/dfi\ > 1. This can 
always be prevented by choosing M large enough. However, from the definition of D, (67), 
we see that D has quite large gradients for r < — y^. We avoid choosing M to be large 
to compensate for this by defining 5^{^) to be (5(/i) and basing the definition of 6~ on 5+ 
with 

6-[^i + 2S+i^i)]=6+i^i)=Sif,), (Cll) 

which is obtained by letting e = 5{iJ, — e) in (C7). In the limit, M ^ oo, (Cll) agrees with 
(C8). 

To summarize: We solve (CI) and hence (66) using the difference scheme, (C3), with 
S+ and 5- defined by (Cll) and (C4). M is chosen to be 25. 
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Fig. 1. Motion of an ion in velocity space, showing the kicks it receives when passing 
through wave-particle resonance. 
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Fig. 2. Comparison of the difference equations, (13), with the Lorentz force law, (2). (a) 
The mapping of the {9, r) plane using (2) with u = 30.23 and a = 2.2. [This is taken from 
Fig. 3(b) of I.] (b) The mapping of the (9, p) plane under T using (13), with 5 = 0.23 and 
A = 0.1424, which is given by (11) with v = 30.23, a = 2.2, and r = 47.5. In each case the 
trajectories of 24 particles are followed for 300 orbits. Thus, in (b) the points, T^{9o, po) for 
< j < 300, are plotted for 24 different initial conditions, {9o,po). 




Fig. 3. Trajectories in the {9, p) plane for small A and (a) 5 = Q, {h) 5 = \, {c) 5 = s/p (not 
equal to or |). The trajectories are obtained by plotting curves for which the Hamiltonian, 
h, is constant. 
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Fig. 4. Eigencurves of the mapping, T, for 5 = and (a) A = 0.3, (b) A = 0.35, (c) 
A = 0.45. (d) The extension of one of the eigencurves in (c) emanating from the hyperbohc 
fixed point with reflection. 




Fig. 5. The iterates of a single point for i5 = and A — 0.3. The point is chosen to start 
on one of the eigencurves emanating from the hyperbolic fixed point and 10 000 iterates are 
shown. Compare with Fig. 4(a). 
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Fig. 7. Plot of h against A for 6 = 0.23, uq = 2, vq = 6. The dashed hne is the asymptotic 
result, (43). 




Fig. 8. The correlation function, Ck, 
Ck was computed using (53) with M 
stochastic region. 



for 8 = 0.23 and (a) A 0.2, (b) A = 0.5, (c) A = 10. 
= 10 and L = 5000. The orbits were all chosen in the 
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Fig. 9. The function, g{A), against A for 5 = 0.11 (crosses), 0.23 (triangles), and 0.47 
(squares), (a) and (b) show two different ranges of A. The point for 5 = 0.47 and A = 0.55 
is at g{A) = 2.8, and so lies off the scale in (a). The function, g, was evaluated using (55) 
and (56) with M = 20, i = 2500, and (a) K' = 100, K = 150, (b) K' = 50, K = 100. The 
solid line gives the approximate form for g, (57). 
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Fig. 10. The correlation function Ck for S = 0.47, A = 0.55, when an accelerator mode is 
present. Here we took M = 100 and L = 5000 in (53). The orbits used to compute Ck all 
lay outside the accelerator mode. 
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Fig. 11. Maximum, minimum, and rms velocities of ions as given by, (a) and (c), solution 
of the exact equations of motion, (2) (taken from I), and, (b) and (d), a Monte Carlo solution 
of (66). In both cases the orbits of TV = 50 particles were followed with a = 20, = 30.23, 
and initial velocity, tq = 23. In (c) and (d) the time scale is altered to show the short time 
behavior more clearly. 
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Fig. 12. The perpendicular velocity distribution function for the particles in Fig. 11 aver- 
aged over orbits 800-1100. (a) The distribution function obtained from the exact equation, 
(2) (taken from I), (b) The distribution function obtained from (66). In each case the 
normalization is such that / 27rr/ dr = 1. 
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Fig. 13. The diffusion coefficient D, (69), as a function of v± for the case discussed in Sec. 
XI. 
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Fig. 14. The steady-state solution for the distributfon function, (a) The two-dimensional 
distribution function, /, as given by the two-dimensional Fokker-Planck equation, (b) The 
one-dimensional distribution function, F(v±), as given by integrating the two-dimensional 
result over (solid line) and solving the one- dimensional Fokker-Planck equation (dashed 
line). 
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Fig. 15. The power dissipated by the wave, Pd, and the power lost to the background 
distributions, Pc, as a function of time for the example given in Sec. XI. The results for the 
two-dimensional and one- dimensional cases are given in (a) and (b) respectively. The units 
of the vertical axes are miV^^noVQ^ ■ The times and Tc are shown in (a). 




Fig. 16. Development of the first-order accelerator mode with m = 1, n = 0, and S 
The cases shown are (a) A = 0.55, (b) A = 0.594, and (c) A = 0.595. For this mode An 
and Au = 0.5927. The crosses show the starting positions of the particles. 




Fig. 17. The range, {Ag,Au), in which first-order accelerator modes exist with (a) s = 0, 
(b) s = 1, and (c) s = 2. 



